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ABSTRACT 

We investigate a massive (S ~ 10 000 gem" 2 at 1 AU) protoplanetary disc model by 
means of 3D radiation magnctohydrodynamics simulations. The vertical structure of 
the disc is determined sclf-consistently by a balance between turbulent heating caused 
by the MRI and radiative cooling. Concerning the vertical structure, two different re- 
gions can be distinguished: A gas-pressure dominated, optically thick midplane region 
where most of the dissipation takes place, and a magnetically dominated, optically 
thin corona which is dominated by strong shocks. At the location of the photosphere, 
the turbulence is supersonic (M ~ 2), which is consistent with previous results ob- 
tained from the fitting of spectra of YSOs. It is known that the turbulent saturation 
level in simulations of MRI-induced turbulence does depend on numerical factors such 
as the numerical resolution and the box size. However, by performing a suite of runs 
at different resolutions (using up to 64x128x512 grid cells) and with varying box sizes 
(with up to 16 pressure scalchcights in the vertical direction), we find that both the 
saturation levels and the heating rates show a clear trend to converge once a sufficient 
resolution in the vertical direction has been achieved. 

Key words: accretion, accretion discs - planetary systems:protoplanctary discs - 
turbulence - instabilities - magnetic fields - MHD radiative transfer 



1 INTRODUCTION & MOTIVATION 

Despite many years of observational and theoretical re- 
search, the physics of protoplanetary accretion discs is 
still not well understood. Significant uncertainties remain 
both with respect to their physical composition as well 
as their tempor al evolution. Their rather short lifetime of 
~ 10 6 to 10 7 yr (|Hartmann et al.lll99Sl ; ISirilia-Aguilar et all 
120061 ) requires a physical process that transfers angular mo- 
mentum outwards efficiently so that the matter can be ac- 
creted fast enough onto the star. Usually it is assumed that 
some sort of turbulence is acting in th e disc, providing an ef- 
fectiv e viscosity that drives accretion (jShakura fc Svunvaevl 
1 19731 . hereafter SS). As of today, the most promising mecha- 
nism for this scenario appears to be the magneto-rotational 
instability (MRI), a linear instability that exists in rotating, 
weakly magneti sed shear flows and operates under very gen - 
eral conditions JBalbus fc Hawlevlll99AlBalbusll2003l . l2009l) . 
Apart from protoplanetary discs, the MRI plays an im- 
portant role also in a variety of other accreting systems, for 
example accretion discs in active galactic nuclei, galactic ac- 
cretion discs and core-collapse supernovae. Although these 
systems cover very different physical regimes, they can be 
modeled using similar numerical methods. The most com- 
mon approach is to employ the shearing-box approximation, 



where one models a small patch of the accreting disc in the 
co-moving frame using special "shear-periodic" boundary 
condi tions which are con sistent with the background shear 
flow (jHawlev et alj|l995l . hereafer HGB). The shearing-box 
model has the advantage that it is easy to set up, and that 
it is not very expensive computationally. In contrast to the 
many works that make use of the shearing-box approxima- 
tion, there are only very few global si mulatio ns of proto- 
planetary discs (IFromang fc Nelson 2009; Flock et aUl20ld : 
iDzvurkevich et alJl201CO . all employing an isothermal equa- 
tion of state. 

The most interesting quantity that numerical simula- 
tions can provide is the magnitude of the turbulent viscos- 
ity, which is usually measured in terms of the alpha pa- 
rameter (SS, HGB). Numerical simulations show that the 
MRI leads to sustained turbulence which transports an- 
gular momentum outwards at a rate that, for the case of 
protoplanetary discs (but may be not for other systems), is 
compatible with observations (|King et al.l 120071 ). The satu- 
ration level of the M RI depends on numer ical parameters 
such as th e box size |johansen et al. 20091). the numerica l 
resolution (Fromang & PaDa loizoull2007l ; ISimon et alj|20091 : 
iDavis et al . 2009; Shi et al. 2010l. heraft er SKH) and the nu- 
merical scheme (jBalsara fc Mever||20101 ). There is also a de- 
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pende nce on the physical p arameters viscosity and resisis- 
tivity (|Fromang et al.ll2007l ) and, to a lesser extent, also on 
the thermodynamics (Flaig et. al. 2009, herafter FK2; see 
also Turner et. al. 2003). When doing MRI simulations, it 
is therefore very important to check the dependence of the 
results on numerical factors and also to include as much of 
the relevant physics as possible. 

Protoplanetary discs are dense systems and therefore 
for a large part optically thick. This means that one should 
include radiation transport in order to properly model the 
transport of radiation energy from the disc midplane, which 
is constantly heated by turbulent dissipation, to the upper 
layers, where the radiation escapes into space. Including 
radiation transport is very important in order to get the 
correct vertical structure, which will be determined by a 
dynamic balance between turbulent heating and radiative 
cooling. It is also an important step towards the goal of 
comparing numerical simulations to actual observations by 
modelling accretion disc spectra, and thereby deriving con- 
straints on physical parameters. 

While such an attempt to include radiation transport 
in MRI-turbulent protoplanetary disc simulations has not 
been reported so far, there exist already several works which 
include radiative transfer in simulations of accretion discs 
that are situated in other phys ical regimes (for example 
iTurner! |2004| ; iKrolik et all l2007t for the case of radiation 
dominated d i scs). T he gas-pressure dominated simulation of 
Irlirose et al.l (|2006l ). hereafter HKS, is most closely related 
to our work, and will be the main source of reference. In such 
simulations it is possible to track in quite some detail the 
path that the energy takes inside accretion discs (HKS) and 
it is p ossible to derive observable consequences (jBlaes et al.l 
l2006t ). 

When dealing with protoplanetary discs one has to 
consider that these are not only dense, but also cold, 
which means that for a large part they will only be poorly 
ionise d, leading to the formation of a "dead zone" (|Gammiel 
Il996l h The effects of a finite resistivity ha ve already 
been studied in isothermal models (for example |Sano et al.l 
|2000| ; ITurner fe Sanoll200Sl ; iDzvurkevich et alj|2010h . In the 
present work, we avoid this additional complication by plac- 
ing our simulation domain near to the star and choosing 
a surface density that is high enough so that the turbulent 
heating provides a sufficient amount of collisional ionisation, 
thereby justifying the assumption of ideal MHD. In a later 
work, we will extend our model by including a physical re- 
sisitvity, thereby arriving at a protoplanetary disc model 
that includes all the relevant physics in a self-consistent 
manner. 

The organisation of our paper is as follows: In §2 we 
describe our numerical methods and the model setup. In 
§3 the results of our simulations are presented. In §4 we 
conclude and provide an outlook on the further research 
that is planned. 



2 MODEL DESCRIPTION 

We apply the stratified shearing-box approximation, i.e. our 
simulations take place in a co-moving rectangular box that 
is small enough so that local Cartesian coordinates (x, y, z) 
can be introduced, which correspond to the radial, azimuthal 



Parameter 


Symbol 


Value 


Radial box size 


L x 


0.08 AU 


Azimuthal box size 


Ly 


0.48 AU 


Vertical box size 


L z 


0.8 AU - 1.28 AU 


Mass of central star 


AL 


1«0 


Distance to central star 


Ro 


1AU 


Surface density 


So 


11356 gem" 2 


Adiabatic index 


7 


1.4 


Mean molecular weight 


M 


2.35 



Table 1. Summary of basic physical parameters of our model. 
The vertical box size of L z = 0.8 AU corresponds to a model 
which covers 10 pressure scale heights in the vertical direction, 
while L z = 1.28 AU corresponds to a model with 16 scale heights. 



and vertical coordinates, respectively. The radial boundary 
conditions are shear-periodic and the azimuthal boundary 
conditions are strictly periodic (HGB). In the vertical direc- 
tion we apply open boundary conditions that allow outflow, 
but no inflow. The simulation box is placed at a distance 
of 1 AU from the central star which has one solar mass. 
We choose a surface density So of order ~ 10 000 gem -2 , 



1700 



g c m given 



which lies in between the value of So 

by the minimum mass solar nebula of lHavashil (|198lT ) and 
the value of So — 50500 gcm~ 2 accor ding to the m ore re- 
cently proposed solar nebula model by iDeschl (|2007l ). With 
this choice of surface density, the temperatures that come 
out of our model are of the order of 1000 — 2000 K in the 
body of the disc, which allows us to describe the the disc 
gas using the ideal MHD equations. The gas is assumed to 
have solar chemical composition, therefore we choose values 
of 7 = 1.4 for the adiabatic index and /j, = 2.35 for the mean 
molecular weight. 

For protoplanetary discs, the stellar irradiation plays an 
important role in determining the degree of ionisation in the 
disc, because stellar X-rays are one of the main ionisation 
sources. However, since for our choice of parameters the disc 
can be considered to be already sufficiently ionised for the 
MRI to work, we may safely neglect this factor. At visible 
wavelengths, the stellar irradiation will only be important 
in the very upper layers of a protoplanetary disc, where ac- 
cording to passively heated models, a temperature inversion 
will occur (|Chiang fc Goldreichl 1 19971 1. In the present work 
we are mainly interested in the temperature profile as it 
results from a balance between active heating by turbulent 
dissipation and cooling by radiation transport. Therefore we 
exclude this region from our simulation domain and neglect 
the stellar irradiation in the visible part of the spectrum 
also. For reference, the basic properties of our model are 
summarised in Table [l] 



2.1 Equations solved 

We calculate the dynamics of the disc gas by solving the 
equations of radiation magnetohydrodynamics in conserva- 
tive fornr. We nrake the one-temperature approximation (i.e. 
we assume thermal equilibrium between matter and radia- 
tion), which means that we need not to solve an additional 
equations for the radiation energy. Radiative diffusion is 
treated within the flux-limited diffusion approach. Includ- 
ing the source terms that arise in shearing box framework, 
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have in our simulations. Four different regions can be distin- 
guished: At the lowest temperatures, in region CD, the dom- 
inant contribution to the opacity is due to ice grains. Since 
these are more efficient scatterers and absorbers at shorter 
wavelengths, the opacity increases with temperature. In re- 
gion ©, ranging from ~ 170 to ~ 200 K, the ice grains melt, 
so the opacity decreases. Region CD is dominated by dust 
grains and the opacity shows only a weak temperature de- 
pendence of the form k oc T ' 5 . Finally, in region ©, the 
temperatures and densities are sufficiently high for the dust 
grains to start melting, leading again to a decrease of the 
opacity. 



Figure 1. Opacity as a function of temperature and density. 
The thick curve gives the opacities for the typical densities and 
temperatures found in our model R7, where M and S correspond 
to the location of the midplane and the surface, respectively, while 
the star * denotes the mean location of the photosphere. For an 
explanation of the different opacity regions © to © see the text. 



the resulting set of equations looks as follows: 
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where the total pressure P to t = p + B 2 
gas pressure and magnetic pressure (neglecting the con- 
tribution arising from radiation pressure), I is the iden- 
tity matrix, / denotes the sum of the gravitational forces 
and the inertial forces arising in the shearing-box system, 
-Etot = vlil - 1) + P v ' 2 / 2 + B 2 /8TC+E laxi is the total energy, 
F = — (Ac/ftp) V-E ra d is the radiative energy flux, and the 
other symbols have their usual meaning. Using linearised ex- 
pressions for the gravitational and inertial forces, / becomes 
(HGB): 



/ = 3pi7 xx — pO z z + 2pQv x z, 



(2) 



with Q the angular orbital frequency. Note that by solving 
the total energy equation, rather than the thermal energy 
equation, all dissipative losses are automatically captured 
and transformed into gas internal energy. In this way the 
heating of the gas is consistent with the dissipation caused 
by the turbulence. 

In the optically thick protoplanetary discs, matter and 
radiation are closely coupled and we can assume that they 
are at the same temperature, which means that the radia- 
tion energy density is given by Erad = fflT 4 , which closes our 
system of equations. For t he flux limiter A we apply t he for - 
mula given in Eq. (28) of iLevermore fc Pomranina (I1981T). 
We us e dust opacities k — n(p, T) as described in lBell fc Linl 
|l99J), which are calculated as a function of the gas density 
and temperature. Fig. [l] shows a plot of the opacity func- 
tion for the typical density and temperature range that we 



2.2 Numerical code 

We use the Cron os magnetohydrodynamics code 
(jKissmann et al.l 120081 ). The code has already been 
successfully applied to simulations of th e turbulent in- 
terstellar medium ( Kissman n et al.l 120081 ). coronal mass 
ejections (|Pomoell et al.ll2008T ) and simulations of MRI tur- 
bulence in accretion discs (FK2). The code solves the ideal 
MHD equation using the conservative scheme described in 
iKurganov et al.l ((200l|). Since we do not include any explicit 
dissipation, the turbulent cascade is cut off at the grid 
scale, where kinetic and magnetic energy are lost due to 
numerical errors and transformed into gas internal energy 
by virtue of the conservative properties of the underlying 
numerical scheme. The fact that the heating is mainly 
due to numerical effects means that there will inevitably 
be some dependence on numerical parameters such as the 
resolution and the details of the numerical scheme involved. 
However, since numerical errors are largest at locations 
where there are sharp gradients in velocity and magnetic 
field, the numerical dissipation should at least qualitatively 
resem ble the actual physical dissipation (HKS, iKrolik et al.l 
120071 ). For high enough resolutions, we expect the results to 
become independent of numerical parameters. It is therefore 
very important to perform runs at different resolutions in 
order to check that the results are actually converged. 

The radiation transport is treated using operator split- 
ting. During the radiation transport step, the equation 
9(eth + E la d)/dt = — V • F is solved, which, using the rela- 
tions eth = phsT/i^f — l)/wrtH and P ra d = aT 4 , can be trans- 
formed into the following diffusion equation for the temper- 
ature: 



dT 
~dt 



(c vP + 4aT 3 )^- = V ■ ^4aT 3 VT, 



Ac 

L 

Kp 



(3) 



where cy = &b/(7 — l)/imH. In order to avoid prohibitively 
short timesteps due to the presence of optically thin regions, 
Eq. ([3j is solved implicitly (for more details on the imple- 
mentation of the radiation transport, see the appendix). 



2.3 Boundary & initial conditions 

In the radial direction, shear-periodic boundary condtions 
are applied, which are consistent with the shear flow in our 
local box (HGB). The azimuthal boundary conditions are 
periodic. In the vertical direction, we use outflow boundary 
conditions which extrapolate the values of the density and 
the velocity in the outermost cell into the ghost zones. In 
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Name 


Resolution 




Box size 


grid cells/Ho 


Orbits 


TD 


(("Rcyn)) 


((«M«w)) 


««» 


16 


32x64x256 


H 


X 6.ffo x 12#o 


21 


200 


isothermal 


0.001 ± 0.001 


0.003 ± 0.004 


0.004 ± 0.005 


15 


32x64x256 


Ho 


x 6H x 10-Ho 


25 


80 


isothermal 


0.003 ± 0.002 


0.011 ±0.006 


0.014 ±0.008 


17 


32x64x448 


H 


x 6.ffo x 14H 


32 


70 


isothermal 


0.004 ±0.001 


0.013 ±0.003 


0.017 ±0.004 


I6D 


64x128x512 


Ho 


x 6H x 12.ffo 


12 


60 


isothermal 


0.003 ± 0.001 


0.012 ±0.005 


0.015 ± 0.006 


R6 


32x64x256 


Ho 


x 6H x 12H 


21 


200 


radiative 


0.001 ±0.001 


0.004 ± 0.003 


0.005 ± 0.004 


R7 


32x64x448 


Ho 


x 6H x UH 


32 


150 


radiative 


0.003 ± 0.002 


0.015 ±0.010 


0.018 ±0.013 


R8 


32x64x512 


Ho 


x 6H x 16H 


32 


90 


radiative 


0.004 ± 0.002 


0.018 ±0.009 


0.022 ±0.011 


R6D 


64x128x512 


Ho 


x 6H x 12H 


12 


60 


radiative 


0.004 ± 0.002 


0.017 ±0.008 


0.021 ±0.010 



Table 2. Overview of the simulation runs that were performed. The box size (third column) is measured in scaleheights Ho refering to 
the initial gas pressure at the midplane. Since the pressure at the midplane does not change very much, this corresponds very well to 
the actual pressure scaleheight during the simulation. The fourth column gives the number of grid cells per scale height in the vertical 
direction. The fifth column tells the number of orbits for which the simulation has been run. "TD" (sixth column) denotes the the type 
of thermodynamics. The last three columns contain the time-averaged values of the turbulent stresses as well as the standard deviation. 



order to prevent inflow of mass, reflective boundary condi- 
tions are applied instead whenever the velocity vector points 
into the computational domain. For the magnetic field, we 
apply vertical field boundary conditions, where B x and B y 
are set to zero in the ghost cells and B z is extrapolated 
from the outermost computational cell. With this set-up we 
do not need any (artificial or physical) resistivity near the 
boundaries in order to prevent inconvenient effects. At the 
boundaries, the temperature is set to a small fixed value of 
Tb = 10 K, which is kept constant throughout the simula- 
tion. This prescription for the temperature allows for free 
emission of the radiation from the disc surface, and it is 
well justified by the fact that in our models the gas at the 
boundary is always optically thin. 

As the initial condition we choose a state of approximate 
vertical hydrostatic equilibrium: 

p(t = 0) = p m exp(-z 2 /2tfo), 

v(t = 0) = Vkcp + 8v, 

T(t = 0) = To, 

B(t = 0) = V x SA, 



(4) 



with the scale height Ho = c g o/fi, where c g o = ^kTo/pmn 
is the initial gas sound speed, and «ke P = — \Qxy is the 
background Keplerian shear flow. For the initial midplane 
density p m and the initial temperature To we choose the 
values p m = 3.93 ■ 10~ gem - (corresponding to a mid- 



plane number density n n 



10 cm ) and temperature 



1500 K, 



respectively. The surface density is then 

2 



To 

So = 11356 gem - ''. In order to avoid excessive shrinking 
or swelling of the disc during the inital phase, the initial 
temperature To has been chosen (by trial and error) such 
that it is roughly the midplane temperature that comes out 
of the simulations for this particular combination of So and 
Ro- Small random perturbations 8v of order 10 -3 c g o are 
added to the velocity. Finally, the initial magnetic field is 
derived from the following vector potential: 

sm(2nx/L 



8A = B - 



2n/L x 



exp(-z /2H )y. 



(5) 



Note that the initial net magnetic flux in the computational 
domain is thus zero, and this property is retained throughout 
the simulation to good accuracy. The value of Bo is taken 
such that it corresponds to a plasma beta of P — 100 at the 
midplane. With this particular choice of magnetic field, the 
whole disc quickly becomes turbulent in less than 10 orbits. 



2.4 Definition of mean quantities 

Since we want to investigate the vertical structure of the 
disc, we are especially interested in horizontally averaged 
quantities. Therefore, for any quantity q, we define the ra- 
dially and azimuthally averaged value (q) xy as 



{q)xy = 



b x L<y 



qdxdy. 



(6) 



In addition to this, we will also use the volume average (q), 
which is defined as 



(«> = 



qdxdy dz. 



(7) 



Lt X LiyLi Z 

The definitions ([6]) and J7| are of course to be understood 
in a discretised sense. 

The most important factor that determines the disc's 
temporal evolution is the amount of angular momentum 
transport that takes place inside the dis c. It is related to 
the r — <fr component of the total stress (jBalbus fc Hawlevl 
Il99a ). which, in the local shearing-box frame, reads 

B x B y 



J- X 



pV X 8Vy 



47t 



(8) 



where Sv y — v y + § Ox. T xy can be decomposed as the sum of 
Reynolds stress Tn, e yn = pv x 8v y and Maxwell stress Tiviaxw = 
-B x B y /4n. 

Following the a prescription of SS, we normalise the 
stresses to the mean gas pressure. We thus define the 
volume-averaged alpha parameter as 

(a) = (aR cyn ) + (aMaxw), (9) 

where 

(cKReyn) = (Tn cyn ) / (P ) and (a M axw) = (Tl&axw}/(P) (10) 

are the volume-averaged Reynolds and Maxwell stresses nor- 
malised by the mean gas pressure. 



3 SIMULATIONS 

Table [2] provides an overview of the simulations that were 
performed. We perform simulations with different box sizes 
and resolutions. The lowest resolution chosen is 32 x 64 x 256 
(Model R6). We also perform a simulation at double the reso- 
lution in all directions (simulation R6D) and two simulations 
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(a) Horizontally averaged azimuthal magnetic field component B v normalised to the initial magnetic field 
strength Bo . The blue curve denotes the location of the photosphere and the green curve is the location at 
which magnetic pressure starts to exceed gas pressure. 




(b) Volume averaged magnetic, turbulent kinetic and thermal energy normalised to the initial gas pressure at 
the midplane. The gas internal energy is about 100 times the magnetic energy (note the different scaling). 




(c) Volume averaged turbulent stresses normalised to the mean gas pressure. 
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(d) Photosphcric temperature at top and bottom boundary. 

Figure 2. Time history of various physical variables for model R7 from to 80 orbits illustrating the changes in turbulent activity and 
the correlations between stresses, magnetic field strength and luminosity of the disc. For further explanation see the text. 



(R7 and R8) with a larger box size in the vertical direction 
and a vertical resolution that is in between model R6 and 
R6D, since recent results indicate that it is the resolution 
in the vertical direction that is critical in determining the 
value of the turbulent saturation level (SKH). In order to 
see how the results from the radiative simulations compare 
to simulations which use an isothermal equation of state, we 
also perform a set of simulations (15, 16, 17 and I6D) which 
have identical physical parameters as the radiative simula- 
tions but which use an isothermal equation of state where 
the temperature is constant kept constant at the initial value 
of T . 



3.1 Time history 

Starting from the initial state described in Sec. l2.3l the MRI 
starts to grow quickly and the non-linear state is reached 
already after a few orbits. After ten orbits the whole disc is 
fully turbulent and remains turbulent for the whole course 
of the simulation. Concerning the time history we now con- 
centrate on model R7. 

The turbulent state is not time-steady, but there are pe- 
riods of high turbulent activity which are followed by periods 
where the disc is less active. This behaviour can be observed 
in Fig. [2j where various physical quantities are plotted as 
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function of time for the first 80 orbits. Fig. 2(a) shows the 



horizontally averaged azimuthal component of the magnetic 
field (which is the dominant component). During an active 
phase the magnetic field is lifted upwards, leading to the typ- 
ical butterfly structure. Recently it has been claimed that 
it is the undulatory Parker instability that drives the up- 
dwelling of the magnetic field (SKH). As one can estimate 
from Fig 



2(a) 



one cycle lasts between 10 and 20 orbits. 
Often the changes in turbulent activity are accompanied by 
magnetic field reversals in one or both sides of the disc. Un- 
like, for example, the solar butterfly diagram, the field rever- 



sals shown in Fig. 2(a) do not follow a regular pattern and 
there is no strict correlation between the magnetic polarities 
in both sides of the disc. 



In Fig. 2(a) we also plot the location of the photosphere 



(defined as the location where the optical depth equals 
unity) and the position of the magnetosphere (defined as the 
location where the magnetic pressure starts to exceed the 
gas pressure). Just between two active phases, the locations 
of the magnetosphere and the photosphere coincide, so the 
magnetically dominated region is optically thin. During an 
active phase, the disc expands due to the stronger magnetic 
forces and the photosphere is pushed outwards. At the same 
time, the magnetosphere is pushed inwards due to the higher 
magnetisation caused by stronger turbulence. Therefore, for 
much of the time, the magnetically dominated region is for 
a large part optically thick. The same phenomenon has also 
been reported in the gas-pressure dominated simulation of 
HKS, although the physical regime is quite different. 

We plot the time history of the energies (thermal energy, 



turbulent kinetic energy and magnetic energy) in Fig. 2(b) 



During active periods the magnetic energy is much larger 
then during a quiet period. The turbulent kinetic energy 
(that is the total kinetic energy minus the kinetic energy 
of the background shear flow) is also larger during active 
phases. In contrast to this, the long-term changes in the 
thermal energy are much smaller than the up-and-down vari- 
ations, which shows that we have indeed reached thermal 
equilibrium. 

The turbulent stresses that determine the angular mo- 



mentum transport in the disc are shown in Fig. 2(c) The 
angular momentum transport is dominated by the Maxwell 
stress which is about four to six times larger than the 
Reynolds stress. This ratio is con sistent with previo us results 
for stratified boxes (for example IStone et al.lll996T ). From a 
comparison between Fig. 2(b) and Fig. 2(c) it is evident that 
the magnetic energy and the alpha parameter are strongly 
correlated. Indeed, it is known that in shearing-box calcu- 
lations, the time-dependence of the alpha parameter can be 
fitted by a formul a of the form (a) ex (B 2 ) + const, (see 
lBrandenburell200Sft . 

Finally, we may ask how the observational appearance 
of the disc will change due to phases of varying turbulent 



activity. To address this question, we plot in Fig. 2(d) the 
photospheric temperature (i.e. the temperature at the loca- 
tion where the optical depth equals unity). If we compare 



Fig. 2(d) with Fig. 2(c) it is obvious that in general a high 



degree of turbulent activity leads to a higher flux of radia- 
tion through the disc's boundaries. The increased flux then 
balances the increased turbulent heating, so that the ther- 
mal energy content in the disc remains nearly constant. This 
observation can be made more quantitative by looking at the 




-5 5 

Lag in orbits 

Figure 3. Cross-correlation between alpha stress and photo- 
spheric temperature for model R7. 



cross-correlation between alpha-parameter and photospheric 
temperature (Fig. [3j : 

C(t) = j[(a)(t - r) - ((a))][T phot (t) - (T phot )] dt. (11) 

Here, (a) denotes the volume- averaged alpha parameter as 
defined in Eq. (|9]) and T p hot is the arithmetic mean of the 
photospheric temperatures at top and bottom. The lag r 
between photospheric temperature and alpha parameter (as 
inferred from the peak in the correlation function) is about 
2-3 orbits. We can compare this to the radiative diffusion 
timescale T ra( j = -H"o/D ra d, where the radiative diffusion co- 
efficient -Dr a d is given by 



D Ta d — 4acT /3cyKp 



(12) 



(cf. the appendix). When taking typical values p — 2.0 • 
KT 9 gcm -3 and T - 1700 K for the density and the tem- 
perature in the region near the midplane (see Fig. [S] and 
Fig. [6] later in the paper), we get r ra d ~ 2.5 orbits, so the 
lag between stress and photospheric temperature agrees well 
with the radiative diffusion timescale. Concerning the tur- 
bulent transport coefficient Dturb = v 2 /fi, we find that near 
the midplane Dturb is typically about one order of magni- 
tude smaller than D ra d , so the energy transport by turbulent 
gas motions is negligible compared to the energy transport 
by radiative diffusion. 



3.2 Vertical structure 

To get an impression of how the turbulent state looks like, 
the reader may first have a look at Fig. |4j where 3D snap- 
shots of some physical quantities are shown for the high- 
resolution model R6D. Panel (a) depicts the magnetic field 
strength measured in Gauss. In panel (b) a plot of the vor- 
ticity measured in units of the angular orbital frequency Q 
is shown. Panel (c) contains a plot of the quantity 



Rz = 



\d z E It 

KpE T ; 



4|&T| 
npT 



(13) 



which is the ratio of the photon mean free path ^ p hot = 1/np 
to the typical length over which the radiation energy (or the 
temperature) varies. 

Concerning the vertical structure, two regions with very 
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(a) Magnetic field strength 
in Gauss. 



\B\ measured (b) Vorticity ui = |V X v\ in units of an- 
gular orbital frequency ft. 



(c) Ratio of photon mean free path to the 
length over which the temperature varies. 



Figure 4. Snapshots from model R6D illustrating the two-layer structure of the disc. The plots in panels (a) and (b) have been taken at 
t = 20 orbits, the plot in panel (c) has been taken at t = 21 orbits. For further comments see the text. (Graphics have been produced 
using the VAPOR visualisation package.) 
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Figure 5. Temperature profiles obtained for the radiative models. 
Averages have been taken from 20 orbits until the end of the 
simulation. 



different physical properties can be disinguished: On the one 
hand, there is the region near the midplane (the region inside 
the first three scaleheights). Here, the magnetic field is ap- 
proximately uniform and of the order of several Gauss. The 
turbulence is subsonic (see Sec. l3.2.3l below) and many inter- 
twined vortices can be observed. Since R z <^ 1 the photons 
are diffusing in this region. On the other hand, there is the 



corona (the region outwards from 4 scaleheights) which ex- 
hibits markedly different features: Here, the magnetic field 
drops sharply, the flow is supersonic and characterised by 
strong shocks, and the photons are to a large part better 
described as free streaming rather than diffusing. 

Due to the properties of the shearing-box setup, on aver- 
age the physical state will appear homogeneous with respect 
to the radial and the azimuthal direction. This suggests to 
consider horizontal averages of physical quantities. However, 
due to the turbulent nature of the flow, the horizontal av- 
erage of some variable at one instant in time shows large 
fluctuations and sometimes strong asymmetries between the 
two halves of the disc. It is therefore necessary to perform 
averages over many orbits in order to achieve a smooth and 
symmetric vertical profile. 



3.2.1 Density & temperature profile 

One of the most interesting quantities that our radiative 
models are able to provide is the self-consistently calculated 
temperature profile. We plot the mean temperature profiles 
for our radiative models in Fig. [5] Except for the low resolu- 
tion model R6 all the radiative models yield nicely consistent 
temperature profiles. In the body of the disc, the tempera- 
ture profile resembles an inverted parabola and becomes flat 
in the optically thin regions in the upper layers of the disc. 
Near the disc midplane, the temperature is about 1800 K. 
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Figure 6. Density profile for models 17 and R7, averaged from 20 
orbits until the end of the simulation. For comparison, the initial 
gaussian density profile has also been plotted (short-dashed line). 
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Figure 7. Gas pressure and magnetic pressure for model R7. Av- 
erages have been taken from 20 orbits until the end of the simu- 
lation. 




The mean position of the photosphere is located at about 
5 scaleheights away from the midplane. The corresponding 
photospheric temperatures are of the order of 500 — 600 K. 
According to the well-known formula 



M = 



8tt aT 4 R% 
3 GM, 







(14) 



this would correspond to an accretion rate M of order 
lO _5 M0yr _1 . The accretion rate can also be estimated 
fr om the a lpha parameter according to the formula given 
in |Pringlel(|l98ll ) 



M = 3nac g HS . 



(15) 



Using c g ~ c g o as well as H » Ho, and plugging in a value 
of a = 0.02 for the alpha parameter (see Sec. 13. 3[) . we get 
again an accretion rate of order lO -5 M0yr -1 , so the two 
results are nicely consistent. 

In Fig. [6] we compare the density profile obtained with 
model R7 to the profile of the isothermal model 17 and also 
to the initial profile. Away from the midplane the density 
profiles flatten due to the additional magnetic support. Since 
for the radiative model the temperature decreases outwards, 
the density in the upper layers is lower as compared to the 
isothermal model. 



3.2.2 Vertical support & Stresses 

In our model, the radiation pressure is small compared to 
the gas pressure almost everywhere. This means that con- 
cerning the vertical support of the disc against gravity, ra- 
diation pressure plays no role, leaving only gas pressure and 
magnetic forces. As can be inferred from Fig. [JJ where gas 
and magnetic pressure are plotted as a function of height, 
the midplane region inside the first three scaleheights is 
gas-pressure dominated. The magnetic pressure is approx- 
imately constant in the midplane region with a slight in- 
crease outwards. Outside the midplane region, the magnetic 
pressure declines exponentially, but not as steep as the gas 
pressure. As a consequence, the disc's corona is magnetically 
dominated. 

The turbulence is not only important with regard to the 
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Figure 8. Turbulent stresses as a function of height for the ra- 
diative models. Averages have been taken from 20 orbits until the 
end of the simulation. 



transport of angular momentum, but also insofar as it con- 
stantly creates new magnetic field and heats the disc gas; 
balancing the losses of thermal energy and magnetic field 
and thus preventing the collapse of the disc. We can deter- 
mine the strength of the turbulence at different locations 
by measuring the turbulent stresses. This is done in Fig. [8] 
where vertical profiles of the alpha parameter are plotted 
for the radiative models. Note that except for the low reso- 
lution model R6 all profiles are very similar, indicating con- 
vergence. We encounter a similar picture as in Fig. [JJ for 
the case of the magnetic pressure: Inside the gas-pressure 
dominated midplane region the stress profiles are roughly 
constant with a slight increase outwards (with the excep- 
tion of the low-resolution model R6, where the stresses drop 
noticably at the midplane). Outside of the midpane region 
the stress profiles decline exponentially. This means that al- 
most all the angular momentum transport happens in the 
midplane region and is there almost independent of height. 
In the corona, angular momentum transport is negligible. 
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Figure 9. Turbulent velocities for model R7, where c g denotes the 
gas sound speed, c g o is the initial gas sound speed at the midplane 
and ctot is the total sound speed (including gas and magnetic 
pressure). Averages have been taken from 20 to 60 orbits. 
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Figure 10. Time-averaged value of the alpha parameter as a 
function of the number of grid cells in the vertical direction per 
scale height for all models. The "error bars" correspond to the 
deviation from the mean value. 



3.2.3 Turbulent velocities 

We now look at the velocity distribution in the disc. In Fig. [9] 
we plot the horizontally averaged turbulent velocity u tU rb = 
v — i>kc P normalised by the gas sound speed as a function of 
height. In the midplane region, the turbulence is subsonic, 
while in the corona it becomes highly supersonic, exceeding 
Mach 5 at the boundaries. For comparison with other works 
we also plot the velocity normalised to the initial isothermal 
sound speed at the midplane and the velocity normalised 
to the total sound speed ctot, where c 2 ot — (p + B 2 /8n)/p. 
Our res ults are in agreement w ith the isothermal simulations 
done bv lMiller fc Stone! (J200CJ) who report Mach numbers of 
about two and with HKS, who report Mach numbers (with 
respect to the total sound speed) between one and two. 

From Fig. [4l middle panel, where a snapshot of the vor- 
ticity is plotted, on can appreciate the highly dynamic and 
tangled nature of the velocity field. In the subsonic mid- 
plane region many intertwined vortices can be distinguished. 
Vortices have been proposed as a means of trapping dust 
particles, thereby helping the formation of larger bodies by 
enhancing the numb er of collisions and slowing down the spi- 
raling into the star (jBarge fc Sommerial 1 19951 ; iTanga et al.l 
1 19961 ). However, the vortices in our simulations are rather 
short-lived, usually lasting much shorter than one orbit, so 
particle trapping in the MRI-generated vortices in the mid- 
plane region will not be efficient. 

In observations, the turbulence will show itself in the 
form of a broadening of spectral lines due to the velocity 
dispersion induced by it. At the mean location of the pho- 
tosphere the turbulent velocities are about Mach 2-3 which 
implies a significant effect on the line widths. The detection 
of CO overtone emission in young stellar objects (YSOs) 
provides a useful diagnostic tool to detect turbulent line 
broadening. The reason for this is that the near overlap 
of CO transitions near the v — 2 — band allows a sep- 
aration of the local broadening (e.g. turbulence) from the 
macrobroadening (caused for example by a disc wind). In 
the case of several YSOs, there is indeed substantial em- 
pirical evidence for supersonic turbulent line broadening of 



the magnitude that we find in our simula tions (jNaiita et al.l 
ll99fj ; ICarr et al.ll2004l ; lHugelmeverll2009l ). 



3.3 Turbulent saturation level 

As has already been remarked in the introduction, in numer- 
ical simulations of MRI tubulence, the turbulent saturation 
level is influenced by a number of numerical factors. By per- 
forming a suite of simulations with different box sizes an at 
different resolutions, we are able to gauge the strength of 
the influence of these numerical parameters. 

In Fig. [TO] we plot the values of the alpha parameter 
as given in Table [T] When increasing the resolution from 
32x64x256 (models 16 and R6), to the double resolution of 
64x128x512 (models I6D and R6D), the values of the stresses 
increase considerably. A similar increase of the stresses is ob- 
served also when only the resolution in the vertical direction 
is increased (models 15, 17, R7 and R8), which is analogous 
to what has been found in the radiative simulations of SKH. 
A look at Fig. [5] suggests that it is especially the region near 
the midplane that is not sufficiently resolved, since the stress 
drops there noticably, while in the better resolved models it 
does not. 

The results depicted in Fig. 1 101 suggest that once a res- 
olution of about 30 grid cells per scaleheight in the vertical 
direction has been achieved, the results will no longer sig- 
nificantly depend on the numerical resolution. All models 
except 16 and R6 are consistent with an a-value of about 
a ~ 0.015 — 0.02, with the isothermal models yielding some- 
what lower stresses that are closer to the lower end of 
a ~ 0.015, while the stresses found in the radiative mod- 
els are grouped around the value a ~ 0.02. No tendency is 
found for the turbulent saturation level to change with re- 
spect to neither the vertical box size nor the resolution in 
the radial and azimuthal directions. 

As a consequence of the lower stress, the heating in 
the low-resolution model R6 is also smaller than compared 
to the other radiative models, leading to significantly lower 
temperatures in this model. Apart from this, the tempera- 
ture profile for model R6 is qualitatively similar to the other 
models. 
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As has already remarked before, the heating does in- 
trinsically depend on numerical effects. However, as can be 
seen from Fig. [5l the temperature profiles for the better re- 
solved models (R7, R8 and R6D) are very similar, suggesting 
that the heating rates are also converged. 
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4 CONCLUSION & OUTLOOK 

We have performed 3D radiative MHD simulations of MRI- 
turbulent protoplanetary discs. Including radiation trans- 
port allows us to obtain a self-consistent picture of the verti- 
cal structure which results from a dynamic balance between 
turbulent heating and radiative cooling. Such simulations do 
thus contain a greater level of realism as compared to pre- 
vious isothermal simulations. Furthermore, radiative simu- 
lations are an important step towards bringing numerical 
simulations into contact with observations. 

In our models two regions can be distinguished: A gas- 
pressure dominated midplane region where the magnetic 
field and stresses are approximately constant on average and 
where almost all the angular momentum transport and dis- 
sipation of turbulent energy take place, and a magnetically 
dominated corona where the stresses and the magnetic field 
are much smaller than in the midplane region. The position 
of the photosphere is highly variable due to changing tur- 
bulent activity. For most of the time, a large part of the 
magnetically dominated part of the disc is optically thick. 
Only at times when the turbulent activity is low do the po- 
sitions of the photosphere and the magnetosphere coincide. 

Since the turbulent saturation level in simulations of 
MRI turbulence may depend on the numerical resolution, 
and since in our simulations the heating of the disc is to a 
large part due to numerical dissipation, it is especially im- 
portant to check what influence numerical parameters such 
as the resolution or the box size have on the outcome of the 
simulation. By running models at different resolutions and 
with different box sizes, we are able to verify previous results 
that it is the resolution in the vertical direction that is crit- 
ical for achieving convergence of the turbulent stresses. We 
find that once a sufficient resolution in the vertical direction 
has been achieved (~ 30 grid cells per scaleheight), neither 
the saturation level nor the heating rates do significantly 
change when increasing the resolution further. No trend is 
found for either of these to quantities to change with the 
vertical box size. 

For our model parameters, where we chose to simulate 
a massive protoplanetary disc, the turbulent heating alone 
already proved to be sufficient to justify the assumption of 
ideal MHD. The logical next step will be to include resis- 
tivity into the model and choose parameters that are ap- 
propriate to a "standard" protoplanetary disc model. When 
including both non-ideal MHD and radiation transport, it 
will be possible to do truly self-consistent simulations that 
do not only include the proper thermodynamics, but also a 
self-consistently evolved dead zone. 
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APPENDIX A: IMPLEMENTATION OF 
RADIATION TRANSPORT 

In this appendix we describe the implementation of radia- 
tion transport into the CRONOS MHD code. By adding the 
thermal energy equation 



de a 

at 



+ V ■ (eth«) 



-pV ■ v - K P p(47tB - c£ rad ) (Af) 



and the radiation energy equation 

^l+V-(E rad w) = -^±V-v+^p(4nB-cE I&d )-V-F, 

(A2) 
we obtain an equation that describes the evolution of the 
combined (gas+radiation) energy: 



d(e th + E Ti 

at 



-+V-[(e th -hE rad ;H = -ptotV-v-V-F, (A3) 



where in general ptot = p+E la d/3. However, since in our sim- 
ulations the radiation pressure is small compared to the gas 
pressure, we neglect it and set p tot — p. Splitting off the radi- 
ation diffusion part an d making use of the flux-limite d dif- 
fusion approximation iJLevermore fc Pomraninelil981i ). the 
radiation transport step becomes: 



d(eth + -Brad) 

at 



-V ■ — V£ rad . 

Kp 



(A4) 



Using the ideal gas law e t h = cvpT and assuming thermal 
equilibrium between matter and radiation, £ rad = aT 4 , we 
close our system of equations and end up with the following 
diffusion equation for the radiation energy: 



, i \ <9-E rad 



iE, 



<)t 



„ Ac„„ 
-V V£ rad . 

Kp 



(A5) 



From this we can estimate the timescale to reach thermal 
balance according to r ~ i? 2 /D rad , where H is the scale 
height and 



D r . 



( e t h 



\4E T , 



+ 1 



Ac 

Kp 



(A6) 



Note that in the radiation pressure dominated case, D rad is 
simply given by D rad = Ac/ Kp. However, for the gas pressure 
dominated case that we are considering, we actually have 
D la d = 4i? rad Ac/ethKp, which means that r is bigger by a 
factor of 4_E rad /e t h and it takes thus much longer to reach 
thermal balance. 

When working in the gas-pressure dominated regime, it 
is natural to replace the radiation energy i5 rad in Eq. ()A5|I 
by aT 4 and then solve the following diffusion equation for 
the temperature: 



(c v p + 4aT' 



1 dt 



V ■ — 4aT 3 VT. 

Kp 



(A7) 
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This equation is straightforwardly discretised in the follow- 
ing manner: 



T, 



«jfc 



T? ]k + 



At 



cvP? jk +MT? jk ) 3 



Ax 2 L i+ij fe *+iife 



d" i ..at: 



1+1 1 



1 r D « 

Ay 2 I % i+h 

Ah[ D[ 



AT" 



ij+3* 



A ^n+1 

iife+i^-'yfc+i 



Z3™. 



,AT" 






where A^. fc = T» + +} fc - T™+\ 
Ac 



i+ijfc 

D" i .. = i 



-4aT c 



K/0 



+ 



i + ljk 



^4aT 3 
up 



ijfc 



(A8) 



and so on. The resulting matrix equation, which involves 
a diagonally dominant matrix, can then be solved by any 
standard matrix solver. For the calculations reported in this 
paper we used the method of successive overrelaxation. 
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